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Abstract 

We investigate longitudinal vibrations of a bar subjected to viscous boundary conditions at each end, and 
an internal damper at an arbitrary point along the bar's length. The system is described by four independent 
parameters and exhibits a variety of behaviors including rigid motion, super stability /instability and zero damp- 
ing. The solution is obtained by applying the Laplace transform to the equation of motion and computing the 
Green's function of the transformed problem. This leads to an unconventional eigenvalue-like problem with the 
spectral variable in the boundary conditions. The eigenmodes of the problem are necessarily complex-valued 
and are not orthogonal in the usual inner product. Nonetheless, in generic cases we obtain an explicit eigen- 
mode expansion for the response of the bar to initial conditions and external force. For some special values of 
parameters the system of eigenmodes may become incomplete, or no non-trivial eigenmodes may exist at all. We 
thoroughly analyze physical and mathematical reasons for this behavior and explicitly identify the correspond- 
ing parameter values. In particular, when no eigenmodes exist, we obtain closed form solutions. Theoretical 
analysis is complemented by numerical simulations, and analytic solutions are compared to computations using 
finite elements. 

Keywords: longitudinal vibrations, viscous boundary conditions, modal decomposition, vibratory response. 



1 



1 Introduction 



In this paper we analyze longitudinal vibrations of a bar with dampers attached at each end as well as at an internal 
point of the bar. This type of problem occurs in modeling structures containing shock absorbers and in control of 
continuous structures with discrete elements. Mathematically, the problem reduces to solving the wave equation 
modified by a Dirac delta term with viscous boundary conditions. When the boundary conditions are classical, 
e.g. the ends are free or clamped, separation of variables or Laplace transforms reduce the situation to a boundary 
eigenvalue problem for a second-order ODE called the Sturm-Liouville problem. These problems are self-adjoint and 
admit a complete system of orthogonal eigenmodes the solution can be expanded into with coefficients determined 
from initial values using the orthogonality. 

When viscous boundaries are present the Laplace transform leads to a boundary value problem with the spectral 
parameter entering boundary conditions. One still gets a system of eigenmodes, but they are not orthogonal in 
the usual inner product, and the eigenvalues are general complex numbers reflecting the non-self-adjoint nature of 
the problem. For some critical values of damping parameters the system of eigenmodes may not be complete, and 
even when it is finding the expansion coefficients in terms of initial data is non-trivial because the eigenmodes are 
not orthogonal. Although studied by mathematicians [TJ [5] such problems and their properties are rather sparsely 
treated in the engineering literature, nevertheless see [3J|31[S1[S] and example 4 in [71 chap 4]. Hull [5] was first to 
treat a bar with a viscous boundary at one end and clamped at the other, but he utilized a non-standard approach 
to decoupling the equations of motions and provided a response only for a harmonic driving force. Udwadia [6] 
appears to be the first to provide a complete closed form solution to this problem via Laplace transform. 

Adding an internal damper at an arbitrary point of the bar, as we do in this paper, significantly complicates 
a closed form solution to the problem. In particular, it is no longer possible to find analytic formulas for the 
eigenvalues since their determination depends on solving algebraic equations of arbitrarily high degree. However, if 
the eigenvalues are found numerically the solution can be written explicitly in a closed form. We obtain the analytical 
solution by taking the Laplace transform and finding the Green's function for the resulting boundary eigenvalue 
problem. In our analysis we are able to take advantage of general mathematical results that simplify calculations 
considerably. For example, although we do find a cumbersome explicit formula for the Green's function it is not 
necessary to find the expansion of its inverse Laplace transform, which depends on eigenmodes and eigenvalues only. 
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Moreover, many qualitative traits of the solution can be gleaned from the characteristic equation for eigenvalues 
directly without computing the vibratory response. Internal damper in a problem with free ends was considered in 

Behavior of the bar is controlled by four dimensionless parameters, the damping coefficients hi and /12 at the 
left and right ends, the internal damping coefficient ^13, and the ratio a/L characterizing the position of the internal 
damper (a is the distance to the damper from the left end of the bar, and L is its length). Since the dimension of 
the parameter space is four it can not be easily visualized. As the parameters are varied, the bar exhibits a variety 
of behaviors including rigid motions, zero damping, super stability and instability. Although a four-dimensional 
diagram can not be drawn we give analytic conditions for all these types of behavior. Much of the unusual behavior is 
due to the fact that we do not restrict /i^'s to positive values they would take if the dampers are realized as dashpots. 
For the negative values we are dealing with so-called active dashpots, or rather 'pushpots', that add energy to the 
bar instead of damping it. Such discrete elements are sometimes used in control problems for continuous structures 

Perhaps the most striking observation is the extreme sensitivity of the eigenvalue distribution to the nature of 
the number a/L. When this number is rational the eigenvalues are generically distributed along p vertical lines in 
the complex plane, where p is the denominator of a/L in lowest terms. This significantly complicates expansion 
into eigenmodes since increasingly larger numbers of them have to be kept. When a/L is irrational this distribution 
appears random. Vibratory response on the other hand, is qualitatively insensitive to the placement of the internal 
damper, and in practice one may want to use ratios with small denominators like 1/2, 1/3, 2/3, etc. A fiip side of 
these observations is that while FEM produces good approximation for the vibratory response at least for short 
times it performs poorly in approximating the eigenvalues. In fact, it produces spurious eigenvalues with large real 
parts that do not converge to actual eigenvalues as the number of elements is increased. When the real parts of the 
spurious eigenvalues are positive FEM will lead to large errors in the vibratory response at large times. 

We organize our presentation as follows. The next section gives the precise problem statement and describes our 
approach towards solving it and the main results obtained. In sections I3I4I we respectively derive analytic formulas 
for the eigenmodes, and reduce computing the eigenvalues to solving an algebraic equation for a/L rational. Section 
[S] discusses in more detail the case, when the damper is placed exactly in the middle of the bar, i.e. a/L = 1/2. 
We compute the eigenvalues explicitly and also give explicit conditions for the undamped behavior of the bar. The 
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Laplace transform of the Green's function of our problem is computed in Section [5] and we discuss its expansion 
into partial fractions. Under generic conditions such expansion exists and is easily Laplace inverted providing 
a convenient way for solving our initial-boundary problem. We also give formulas for special combinations of 
parameters when the Green's function can be inverted analytically. Section [7] presents theoretical analysis of 
eigenmode completeness for our problem, and discusses the physical meaning of critical cases when this completeness 
is lost. In section |5] we use the eigenmode expansion of the Green's function to write the vibratory response of the 
bar to initial data and external force. Section |5] compares analytic and FEM solutions in several parameter regimes 
and discusses spurious eigenvalues produced by FEM. Finally, we draw our conclusions in sectionlHl Appendix gives 
derivations of some formulas used in the main text. 

2 Main results 

We begin with the problem statement. Figure [1] depicts a bar of length L free to move horizontally suspended by 
two dampers at each end and by one at the distance a from its left end. Symbols p, A and E represent the density 
of the bar, the constant cross-sectional area and its modulus of elasticity respectively, the wave speed along the 
bar is denoted c := {E / p)^/"^ . Let ci, C2 and C3 be the damping coefficients of the left, right and internal dampers 
respectively, we set hi :— -^ci, := -^02 and /13 := jg^cs (the extra 1/2 simplifies some formulas). These hi 
along with a/L are the dimensionless parameters that determine qualitative behavior of the bar. Since in our case 
the bar can move rigidly just as in the problem with free ends, we write the equation of motion in the absolute 
frame that remains at rest at all times. At t = the left end of the bar is assumed to be at the origin, and u{x, t) 
denotes the displacement of the point with initial coordinate x at time i, see FiglTJ 



u{x, t) 




Figure 1: A bar with viscous ends and internal damper. 
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The system is described by a modified wave equation 



utt{x,t) + 2/i3C(5(a; - a)ut{x,t) - c^Uxx{x,t) = p{x,t), (1) 

with the boundary conditions 

UxiO,t)-^utiO,t) = and u4L,t) + ^ ut{L,t) = 0. (2) 

Here p{x^ t) is the external force per unit mass, and the subscripts x, t denote partial derivatives with respect to space 
and time. Given u{x^ t) the solution in the frame that moves along with the left end of the bar is u{x, t) — u(0, t). 

To solve the problem we use Laplace transforms. Setting the external force and initial data to zero and taking 
Laplace transform we get a homogeneous boundary problem for U{x, s) :— C[u{x^t)\ 

2 

^ U{x, s) + 2/i3- 5{x - a) U{x, s) - U^xix, s) = 0, (3) 



UJO,s)-hi-U{0,s) = and Ux{L, s) + h^- U{L, s) = 0. (4) 
c c 

We will explicitly compute the Green's function G{x,^, s) for this problem and use it to solve the original initial- 
boundary problem. 

Unfortunately, the inverse Laplace transform of G{x, f , s) can not be computed in closed form except in special 
cases. To invert the Laplace transform for G{x, ^, s) we use the spectral method. Note that one gets the same 
boundary problem by separating variables and looking for solutions of the form u{x,t) ~ ipa{x, s)e'^* , where s is 
the spectral parameter. System is almost an ordinary boundary eigenvalue problem except for the presence 

of s in the boundary conditions. One expects that it is solvable only for special values s„, the eigenvalues, with 
faix, Sn) being the vibrational eigenmodes of the bar. This is indeed the case, and for a/L — q/p (with integers q,p 
in lowest terms) they can be grouped into p series sit\ one for each root of a degree p algebraic equation. These 
roots are the only quantities to be computed numerically. Once they are determined, si*"' and the eigenmodes can 
be written explicitly. 
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We then show that under generic conditions the Green's function [G(x, ^, s)] can be expanded into a series 
over the eigenmodes and give an expUcit formula for the expansion coefBcients Aii \ With all the pieces in place 
we derive our main result, a series solution for the vibratory response of the bar with initial data u{x,0), u(x,0) 
subjected to the external force p{x,t): 



u{x, t) 



1 



hi + h2 + 2h3 



fc=l ri= — oo 



hiu{0,0) + h2u{L,0) + 2h3uia,0) + - / ?i(^,0)+ / p(f , r) dr 

cJo L Jo -I 

c hi u{0, 0) + ch2 u{L, 0) ipaiL, 4''^) + 2c /13 u{a, 0) ipaia, s^f^ 



(5) 



Jo 



The other major result is a complete analysis of critical behavior. Namely, we prove that the eigenmodes and 
associated modes are sufficient to expand the Green's function if and only ii hi 7^ ±1. When /12 = 1, /13 = or 
hi — h2 ^ 1 we derive alternative solutions in closed form. 



3 Eigenmodes 

In this section we compute the eigenmodes of our problem analytically. To aid notation and understanding we 
first compute the eigenmodes for the simpler system without the internal damper (^13 — 0). Define (p{x, s), ip{x, s) 
as solutions to ^U{x,s) — Uxx{x,s) = satisfying only the left and the right boundary condition from Eq.(l4]) 
respectively. This defines them up to a constant multiple and we make them unique by normalizing if{0, s) = 1 = 
tp{L,s). One easily finds 

ip{x, s) = cosh(^ — j + hi sinh(^ — j = 2 [(^ + hi)e^ + (1 - hi)e~^] (6) 



■>p{x, s) — cosh^^^^^ — + /i2 sinh^ 



' s{L — x) 



s(I.-x) 3(L-x) 

(1 + /i2)e^^ + (1 - h2)e — 



(7) 



By construction, any solution to j^U — Uxx = satisfying the left (right) boundary condition must be a multiple 
of ip{x, s) ( V'(a;, s) ). Since the eigenmodes are supposed to satisfy both we conclude that they only exist for values 
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of s that make ip{x, s) and ipi^i ■'^) proportional. For such values either one of them can be taken as the eigenmode. 



Proportionality happens if and only if the Wronskian 
eigenvalues — 0. For future reference we denote 



'fix ^'x 



is zero giving us an equation for 



A(s) := -"^Wlip, V] = (1 + hih2) sinh(^ 



{hi + h2) cosh{^^ 



1 r 

2 



(1 + hi){l + /i2)e^ - (1 - hi){l - h2)e~ 



(8) 



and remark that the eigenvalues other than 5 = satisfy A(s) = or 



2^ {l-hi){l-h2) 

e " = 



{l + hi){l + h2) 



The latter can be solved explicitly: 



c 

2L 



In 



{l-hi){l-h2) 



{l + hi){l+h2) 



+ i 



Arg 



(l-/ll)(l~/l2) 
(l+/ll)(l+/l2) 



2?™ 



(9) 



Substituting s„ into Eq.® gives the eigenmodes (/5(a;,s„), also explicitly. 

We now apply the same scheme to Eq.(l3]) when hj, ^ 0. As above, denote by ipa{x,s) {ipa(x,s)) solutions 
to Eq.Q satisfying the left (right) boundary condition and normalized to be 1 at the corresponding boundary. 
Consider ipa first. Since the damper at x — a does not affect the equation on [0,a) we have ipa{x,s) — ip(x,s) on 
this interval. 

2 

For X > a our (pa again satisfies U ^ Uxx = but there must be a jump in its first derivative at a to 
produce the 2/13^ — a)U term in Eq.(|3]). Along with continuity at a we have ipa{a,s) — (p{a,s) — and 
(/?J,(a, s) — ip'{a, s) — 2/13^93(0, s). Since the difference ipa^ (p also satisfies — Uxx = we see by inspection that 
(/3a (a;, s) — ip{x, s) + 2/13 p{a, s) sinh ^ ''(^^ °) ^ for x > a. 

One can compute ipa analogously or notice that by symmetry it can be obtained from ipa by changing x to L — x, 

{1, .T > 
the X < a and 
0, a; < 
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X > a cases can be unified into 



j 



(10) 



i^aix, s) = ip{x, s) + 2/i3 H{a — x) "ipia, s) sinh^- 
To further aid in notation we define 



(11) 



Aa{s) ^ --W[ipa, V^a] = A(s) + 2/i3<p(a, s)V'(a, s) 



(12) 



and compute explicitly 



Aa(s) = {1 + hih2 + h3{hi + h2)) smh(^^^ + (hi + h2 + hsil + hih2)) cosh(^^^ 
+h3{h2 - /ii) sinh( ^^^~"h + hsil - hih2) coshf " ''^ 



\ c 



\ c 



(13) 



We will mostly use the exponential form of this expression 



Aa{s) = i [(l + /^l)(l + /^2)(l + /^3)e'^^-(l-/^l)(l-/^2)(l-/^3)e- 
+(l - hi){l + h2) h3 e^^^^ + (1 + - /12) hs e"^^^^ 



(14) 



Eigenvalues are the solutions to sAa{s) = 0, but Aa(s) contains four distinct exponents and Aa(s) = in general can 
not be solved explicitly. Still, if the zeros are found numerically one can get explicit expressions for the eigenmodes 
from Eq. pl]| . Note that s = corresponds to a rigid displacement of the bar since (pa{x,0) = 1. We address 
computation of other eigenvalues in the next section. 



4 Eigenvalues 

We now need to compute the eigenvalues, i.e. solve the characteristic equation sAa(s) — 0. As we saw in the 
previous section, when there is no internal damper the answer is given explicitly by Eq.®. This is possible because 
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only two different exponents appear in A(s). However, in general Aa(s) contains four different exponents and there 
is no analytic formula for its zeros. Nonetheless, it is possible to reduce this transcendental equation to solving an 
algebraic one. Let us multiply the exponential form of by and set: 

Ai^il + hi){l + h2){l + h3), A2^h3il~hi){l + h2), A3 =/l3(l + /^l)(l-/i2), A4 = - ( 1 - /ll ) (1 - /la) (1 - M • 



In terms of a new variable ^ :— 2sL/c the equation becomes 

Aie^ + Aae^^'*)^ + Aae*^ + ^4 = 0. (15) 



The left hand side is a so-called exponential sum with real exponents, its zeros can be distributed in complicated 
ways. However, these complications only arise when the exponents in Eq. ljlSp are incommensurable, i.e. a/L is 
irrational '9, chap VI]. But irrational numbers can be approximated by rational numbers with arbitrary precision, 
so one can assume, for all practical purposes, that a/L is rational. We do so from now on. 

£ 2sL __ 

Let a/L = q/p in lowest terms with positive integers p,q. A further substitution z — ep = e~ reduces Eq. (jl5p 
to an algebraic equation 

AizP + A2zP-'^ + Asz"^ + = 0. (16) 

li Ai ^ it has p roots zi, 2:2, Zp counting multiplicity, and we need not deal with infinitely many zeros of 
Ea. ((T5|) directly. If moreover A4 ^ none of these roots is 0. When Ai or A4 do vanish we get critical cases that 
require special treatment. 

In general, roots Zk can not be found analytically as functions of hi. However, there is a standard way of setting 
up a matrix with the characteristic equation Eq. (|16p so that they can be found numerically as its eigenvalues using 
MATLAB or Maple. Once Zk are found, the eigenvalues can be expressed as (cf. Eq.®) 

4f' = |^[lnkfc|+z(Arg(zfe) + 27rn)], n G Z, fc = 1, . . . ,p . (17) 

Figl2] shows a typical distribution of eigenvalues, MATLAB function eig() was used to compute the roots z^. 

If all Zk are distinct then si*^"* fall into p infinite sequences with Re(s) = f§ In \zk \ equispaced along vertical lines 



9 



s 
« 

I -2.5 



» 

-4-5 



80 p 

60 - 

40 - 

20 - 

— 

-0.5 

-20 - 

-40 - 
-60 - 
-80 - 



-3.5 




Figure 2: Distribution of eigenvalues for a/L — 2/5 (left), and a/L — 41/100 (right) with hi — 0.3, /12 = 0.9, 
/13 = 0.7, c = 0.3 and L = 1.8. Horizontal lines are the real axes, vertical lines are parallel to the imaginary axes. 

in the complex plane. The lines are not necessarily distinct even if the roots are, since the latter may have equal 
absolute values. In fact, if Eq. (|16p has a complex root Zk then its conjugate Zk is also a root, and the corresponding 
line contains two different with arguments of opposite signs. The eigenvalues are spaced twice as 

densely along such lines, see Figl2] If more roots have the same absolute value the density increases further, but 
the eigenvalues remain distinct as long as the roots are. 

However, if Eq. (jl6l) does have multiple roots then according to Ea. (jl7p we get an infinite sequence of eigenvalues 
with the same multiplicity. Multiple eigenvalues create extra difficulties in computing the Green's function since 
the eigenmodes are no longer sufficient to express the response (the associated modes also enter). We shall assume 
henceforth that all the roots of Eq. p^ are simple. Then all the eigenvalues are also simple with the possible 
exception of s = 0, which will have multiplicity two if Aa(0) ^ hi + h2 + 2 ^ 0. It corresponds to a bar moving 
as a rigid body with constant velocity. 

When p is small Zk can be found in closed form. The simplest such case is p = 2, g = 1 so that a = ^, i.e. the 
internal damper is placed exactly in the middle of the interval [0,i]. We wish to analyze this case in more detail 
next. 
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5 Damper in the middle 

When a/L = 1/2 Eq. p^ reduces to a quadratic equation: 



Az'^ + Bz + C ^0 with (18) 

A={l + h,){l + h2){l + h;), B = 2h3{l-h,h2), C=-(l-h,){l-h2){l-h3). 

As above, we exclude critical cases so that A and C 7^ 0. Then Ea. ([T5| has two real roots, two complex 
conjugate roots or one real double root according to whether 

D ■■= = (1 - hl){l ~ hi) + hi {K - h,f (19) 

is positive, negative or zero respectively. Moreover, we have explicitly 

^ -h^{l-h,h2)±^ 

^'^ {l + h^){l + h2){l + h^)■ ^ ' 

Accordingly, for D > we have two sequences of eigenvalues spaced along two distinct vertical lines, see Ea. (IT7|) . 
They merge into a single line of double eigenvalues for 13 = 0, and for D < we have both sequences half-spaced 
along a single line with 



Re(s) = ^ln|^i,2| = ^ln 



(l-/ii)(l-/i2)(l-/»3) 



[l + hi){l + h2){l + hs) ■ 

Even if D 7^ when /ii + /i2 + 2 /13 = we have a double eigenvalue at zero because of the rigid mode, in which 
case D^{hl + hl- 2f ^ 0. 

As an application, we will determine conditions for 'undamped' behavior of the bar. The idea is this: when 
hi 's are negative they describe active dashpots that add energy to the bar instead of draining it [6 . It is therefore 
possible that for some combinations of values the same amount of energy is being added as is being drained by the 
dashpots - the bar behaves as if it were undamped altogether. Clearly, no damping means that all the eigenvalues 
are purely imaginary or equivalently, that all the roots of Ea. (jl8|) lie on the unit circle since z = . 

If both roots are real then z = 1, 1, z = —1, —1 or z = 1, —1. This leads to C ~ A, B = ±2A or C = —A, B = 
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respectively. If one of them is complex then its conjugate is also a root and z = e^'^ with ^ ^ 0, tt. We have 

Az^ ^Bz^C = A{z - e^){z - e'") = Az^ - 2Acos6' + A, 

so C = j4 and cos^ = —B/2A. The latter condition can be satisfied by a ^ ^ 0, tt if and only if \B\ < 2\A\. Two of 
the real cases above will be subsumed here if we allow \B\ < 2\A\. 

To make these conditions explicit in terms of hi's let us start with the case /13 = 0. This automatically means 
i? = 0, in particular |i?|<2|A|,soC = A and C = — A give us the two available possibilities: 



(1) 



/13 = Us = 

and (2) <^ 

/I1/I2 = -1 /li+/l2 = 0. 



In each case hi or /i2 can be chosen arbitrarily and the other parameter is uniquely defined. When /13 7^ we have 
a similar situation for C = —A: 

hih2 = 1. 



(3) 



Note that in the second and third case we have hi + h2 + 2/13 = 0, which means that there is a double pole at 
s = 0. Therefore, we have not just zero damping but rigid motion: the bar will be moving with constant velocity 
in addition to oscillations. 

When ha ^ the case of C = A becomes unexpectedly complicated. After cancellations C = A reduces to 
hih2 + hihs + /12/13 + 1 = and we also get 

1 - hlhl 

cosO — 



{l-hl)il-hl) 

We can solve for hs since /ii + /12 = leads to one of the above cases, but writing out \B\ < 2\A\ is not very helpful 
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for determining /ii,/i2- We will assume instead that one picks hi and cos 6* and then solves for h2,hz: 



(4) 



hi+h2 



, I /(l+COSe)-/l{ COS0 



One can see that the range of parameters in case (4) is two-dimensional in contrast to the first three cases. Still, 
the choice of hi and cos 6* is not entirely arbitrary. Aside from the obvious restriction — 1 < cos 9 < 1 one has to 



pick hi so that the expression under the square root is non-negative. For example, if cos = 1 then /12 — ±^^2 — hf 
and we can only pick hi that satisfies \hi \ < V2- In particular, \hi\ < 1 always works for cos9 > 0, and \hi\ > 1 
for cos 9 < 0. 

For C = A and B = zL2A, which corresponds to cos 6* = ±1, we have multiplicity in eigenvalues. Indeed, Eq. p7| 
shows that all eigenvalues on the imaginary axis are now double zeros of Aa(s). Therefore, we have oscillations 
with linearly increasing amplitudes. Moreover, when cos 9 = 1 we have that is even a triple pole of sAa(s) and 
hence a triple eigenvalue. Physically, this means that the bar does not just move rigidly, it is accelerating. 

6 Green's function 

By definition, the Green's function G(x,^, s) for system ©-dH) satisfies 

2 

^ G + 2/i3- 6{x - a) G - G,,(x, s) = <5(x - C), (21) 



G,(0,^,s)-/ii-G(0,f,s) = and G,{L,^,s) + h^- G{L,^, s) ^ 0. (22) 
c c 

It can be computed along the same lines as ipa, ipa in section |3l and has different analytic expressions depending 
on relative positions of x, a and ^. Consider the case a < ^ first. As a function of x, G satisfies Eq.Q for a; < <^ 
and X > ^. Therefore, it is equal to Aipa{x, s) on [0, ^) and Bip{x, s) on (^, L\. At x = ^ it is continuous, but has a 
jump in the first derivative to produce 5{x — ^) in Eq. pTj) . Namely, G2;(^+, ^, s) — Gx{^^ s) — —1 because Gxx 



13 



enters Eq.dH]) with minus. Therefore, A, B can be found from the system 



Solving for them in the matrix form we get 



-V'a fa] 



But for ^ > a we have from Eq. pTj) that i^aii, s) = s) so that W[(pa, ipa] = M^[</5a, V"] 
Therefore, 



':Aa{s) from Eq.dl 



sAa{s) 



and B = c 



SO that 



sAa(s) 



<y9a(a;,s)?A(^,s), X <^ 



(23) 



Anafogously, for a > ^ in the latter case we have 



sAais) 



(24) 



It will be convenient for us to rewrite G in a form that is both more explicit, and makes its symmetry G(a;, ^, s) 
G{^,x,s) manifest. To this end, we introduce 



ip{x,s)'ip{^,s), X < C 
ip{^,s)ip{x,s), x> ^ 
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and compute 



9v^{x, £,,s) = - [(1 + hi){l + /i2)e = + (1 - /ii)(l + /i2)e ^ 



+(1 + /ii)(l - /i2)e 5 + (1 - /ii)(l - /i2)e ^ 



(25) 



Analogously, let 



sinh ^(a;, s), a; > ^ 



and gsip{x,£,,s) := < 



sinh 2^2^ ip{x, s), X <^ 



Then 



9si>{x,^, s) 



s(L-a-lx-il) s(L + a-(x+0) 

(1 + h2)e 5 - (1 + h2)e s 



+ (1 - /i2)e (1 - /i2)e s (26) 



^^^(aJ.CjS) = 



(l + /ii)e s -(l + /ii)e s +(l-/ii)e s - (1 - /ii)e ^ 



Since the g^^ part is common to all arrangements of a;, a and ^ we get 



• (27) 



G{x,^,s) = 



sAa(s) 



9<pip {x, ^, s) + 2h3 H{x - a)H{^ - a) ip{a, s) Qs-^, {x, s) 
+ 2/i3 H{a - x)H{a - ^"(0, s) gs,pix, ^, s) . 



(28) 



Note that the last two terms are non-zero only when x and ^ are on the same side of a. Therefore, whenever a 
separates x and ^ the Green's function reduces to the first term. 

As mentioned earlier, to solve the original problem we need to invert the inverse Laplace transform T{x, ^, t) := 
£~^[G{x,^,s)], but G is too complicated to allow inversion in a closed form. We are forced to expand it into a 
series over functions with simpler dependence on s and invert it termwise. In the spectral method, which we follow 
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here, these functions are the partial fractions j^^Jpjm and C^^ (s-p)^ ~ {m-iy. ''^here p are the poles of G. 
Since the denominator of G is sAa{s) its poles are exactly the eigenvalues si'^'' from Eq. (fT7|) and (none of them 
cancel with the numerator). 

Let us assume for the moment that expansion into partial fractions is possible and moreover, that all the poles 
except possibly s — are simple. A general theorem on Green's functions (TUJ chap 1, sect 3] implies that 

G(a;,^,s) = Go(x,e,.s)+ }^ ^y- , (29) 

sL'^Vo yS — Sn ) 

where Gq is the principal part at s = 0, (faix, sl'^^') are the eigenmodes corresponding to si^\ and A^'^^ are numerical 
coefficients. As shown in the Appendix, in our case the latter can be computed explicitly: 



A^' = ''-^l ^je,sW)^d^ + %.(0,sW)2 + %.(i,sW)2 + ^^.(a,s«)^ (30) 



Jo c c c 

If s = is also a simple pole, i.e. Aa(0) 7^ 0, then Go(x,^, s) = ci^(x, ^, 0)/sAa(0), where F is the bracketed 
expression in Eq. ipS)) . By inspection, from Eqs.(|5)).([7)). piI]) .([TT |) we see that = = = V-'a = 1 for s = 0. Recall 
also that Aq(0) = hi + /i2 + 2/13. Therefore, F(a;,^,0) = 1 and 

The Laplace transform is now easily inverted and 

p 00 

nx,S„t):=L-'[G{x,ts)]= , ,^2/. + E E -(^ ^a(x, s^) V'a(C, ^i"^) e^"''* (32) 

If s = is a double eigenvalue, i.e. a simple zero of Aq, the answer is more cumbersome. For example, when 
hs = —hi^li. and /12 = 7^ all the eigenvalues are on the imaginary axis and there is a double pole at zero. The 
principal part of the Green's function at s = is provided with Gq{x, ^, s) = — + % where Ciand C2 are computed 
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in the usual way [TTl chap 5, sect 5-11]. If, for convenience, we set 



1 , a is on the same side of x and ^ 



Ha{x, e) := H{x - a)H{^ ~ a) + H{a - x)H{a - 



(33) 



0, a separates x and ^ 



then ciand C2 are provided as 



hic 



{hl + \){\x-^\-{x + + '^a)Ha{x,£,) 



(34) 



ci = 



2{L{hl-l)-a{h\-l)) 



{hi + l)|a; - ^1 + {hi - l){x + + 2L 




(35) 



C2 = 



L{hl-l)-a{h\-\) ■ 



The corresponding terms in V{x,(^,t) are ci + C2t. Physicahy, the bar is moving as a rigid body with constant 
velocity in addition to vibrating. When s = has higher multiplicity the bar will be accelerating. In the case when 
infinitely many eigenvalues are multiple, the simple template Eg. (1291) for the partial fraction expansion no longer 
applies. We investigate when such expansion is possible at all next. 

7 Critical cases 

Through the previous section we assumed that the Green's function can be expanded into partial fractions over its 
poles. This is not always the case as observed already in [6^. In general, if one takes a function like e^'* with no 
poles at all, it will have no expansion in terms of partial fractions. The usual way of proving convergence is to apply 
the Cauchy residue theorem to circles of increasing radii. But for this argument to work, the contour integrals over 
the circles must tend to zero as the radii tend to infinity. This is violated for functions like . To ensure the 
desired convergence it suffices, for example, that |G'(a;,^, s)| < const/|s| outside disks of any fixed size surrounding 
the poles of G (the constant will depend on the size chosen). When this inequality does not hold we call the case 
critical. 

According to Eq. ipS]) . G{x^ ^, s) = ^"^^^^j^ with F{s) being the expression in brackets. Since 1/s factor is already 
present it would suffice that F{s) / /S.a{s) be bounded away from the poles. A quick look at explicit expressions for 
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F and shows both to be exponential sums of the type aie"^^ + • • • + ams"'"^^ with real exponents a^. Quite a 
bit is known about such sums, e.g. their zeros are located within a vertical strip |Re(s)| < w, and within this strip 
the sum is bounded. Moreover, on the complement to all disks of any fixed size surrounding the zeros exponential 
sums are uniformly separated from Consequently, we only need to worry about F{s)/ Aa{s) being bounded 
when Re(s) — > ±cx). 

Clearly, for Re(s) — ^ oo (—00) the term with the largest (smallest) ai dominates an exponential sum. We 
conclude that F/Aa is uniformly bounded away from the poles if and only if the largest (smallest) exponent in 
the denominator is greater (less) than the largest (smallest) exponent in the numerator. To put it differently: for 
boundedness all exponents in the numerator must lie between the largest and the smallest ones in the denominator. 
Let us take count of these exponents in Eg. ([25)1 for £, > a, the other case is analogous. We find that in the numerator 
the following terms occur 

ef(i-l-"«l), et(^-(^+«), e-*(^-l^-«l), e-f(^-(^+«» (36) 

g|(L-2a-|2;-4|)^ (L+2a-(x+4 )) ^ f (L-2a- Ix^CI ) ^ ^-j.{L+2a-{x+i)) ^ ^3^^ 

and the terms in the second row are multiplied by when x < a. On the other hand, Aa contains 

We see that the exponents of Eq. p6|) do lie between —L and L, so the condition for boundedness is satisfied assuming 
that coefficients in front of e^^ and e^i^ in Aa{s) are non-zero. By Eq. p4)) this means that (l + /ii)(l+/i2)(l+/i3) 7^ 
and (1 — /ii)(l — /i2)(l — /13) 0, or equivalently hi 7^ ±1 for i — 1,2,3. For ODE of any order with linear boundary 
conditions for non-criticality are derived in [2]. 

When say /12 is 1, the only surviving exponents in Eg. (1551) are the first and the third. This changes the calculation 
since some of the exponents in the numerator, like L — {x + are perfectly capable of being less than L — 2a. One 
may hope that 'bad' terms in the numerator also vanish, but no, for /i2 = 1 the first two exponents from each row 
of Eo. ((55|) have non-zero coefficients. 

When hi = ±1 for at least one i the spectral method no longer applies and one has to look for other ways to invert 
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the Laplace transform. Before discussing them let us explain physical reasons for critical behavior. Take h2 — I and 
note that in the original initial-boundary problem the right boundary condition becomes ut{L,t) + cUx{L,t) = 0, 
while the equation away from the internal damper is utt — (? Uxx — 0. It is well-known that its solution splits into 
left and right traveling waves which satisfy Ut — cux = and ut + cux = respectively. The left traveling waves play 
no role near the right boundary (there is nowhere for them to come from), while the right traveling waves satisfy 
the boundary condition automatically. In other words, the right boundary condition has no effect whatsoever, and 
we effectively get a problem on the semi- infinite interval [0, cxd), but with initial data restricted to [0, i]. Physically, 
the right boundary becomes transparent allowing the waves to pass through it without any reflection back. If in 
addition there is no internal damper {h^ = 0) then standing waves, a.k.a. eigenmodes, can not form, so of course 
one can not expand over them. When the internal damper is present, it reflects some of the waves at a forming 
standing waves on [0, a], but not on [a, L]. The expansion is still impossible even though some eigenmodes are 
present. More physical analysis of transparent boundaries is given in [Sj. 

When /i2 = — 1 the situation can be analyzed similarly. Now the boundary condition is ut{L,t) — cUx{L,t) — 
while the equation prescribes Ut -f cux = 0. This is only possible if Ut{L,t) — Ux{L,t) — 0, conditions impossible 
to satisfy in general. If the initial disturbance is localized away from the right boundary the solution will exist only 
up to a time necessary for the right traveling waves to reach it. As ft,2 approaches —1 the reflection coefficient at 
the right boundary grows without bound [5], so the physical reason for the non-existence is that reflection would 
involve an infinite energy transfer. One can not say however, that solution blows up in finite time since it behaves 
'normally' before the boundary is reached and ceases to exist at that instant. 

How does one invert the Laplace transform in critical cases? Suppose first that there is no internal damper 
(/13 — 0) and the right boundary is transparent (/12 = 1). Then there is only one exponent left in the denominator 
of the Green's function Aa(s) = (1 -I- hi)e~ . This exponent can be combined with the ones in the numerator so 
that the Green's function reduces to a linear combination of exponents divided by s 



Since C ^[■^e °"'] — H{t~ a), where H is the Heaviside function, r(a;,^,t) :— C ^[G'(x,^, s)] can even be found in 



G(x,e,s) 




1 - hi 

l + hi 



(39) 
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closed form: 



r(x,c,i) 



H{ct -\x- ei) + ^-5^ H{ct -{x + 0) 



(40) 



1 + hi 

Physically, we get left and right traveling waves with the former reflected from the left boundary just once. Thus, 
it is not surprising that the solution is a finite superposition of traveling waves, just as it is for the wave equation 
on the entire line according to the d'Alembert's solution. 

Summarizing, for ~ we have a simple dichotomy: either both boundaries arc non-transparent and the 
solution can be found by the spectral method, or one or both are and the solution can be found in a closed form as 
a finite superposition of traveling waves. When ^ this dichotomy fails because standing waves may be able to 
form on part of the interval between one of the boundaries and the internal damper. Then one is forced to either 
combine traveling and standing waves, or to use an infinite number of traveling waves. We shall not treat such 
intermediate cases in this paper. 

However, if both boundaries are transparent (hi — h2 ^ 1) the standing waves can not form at all and one can 
find a closed form solution again. Then from Eq. 



2(l + /i3)s 



2(l + /l3) 



Hid - |x - CI) + hMx, (H{ct -\x- i\) - H{ct -\x + i-2a\ 



(41) 

(42) 



8 Vibratory response 

When computing the Green's function we set all the initial-boundary data to zero and it is now time to bring it 
back. The Laplace transform of the original system ([T]),® is 



^ TK \ , Oh ^ Xf \TT( \ TT ( \ U{X,Q) U{X,Q) +p{x,s) ^h'i 

-ir (7(x, s) -f 2/13- d(x - a) U(x, s) - Uxx[x, s) = s 5 1 ^ h 2 — d{x - a) u[x, 0), (43) 

c c c c c 



U40,s)-hi-U{0,s) = -— u(0,0) 
c c 



and UJL,s) + h2-U(L,s)^—u(L,0). 

c c 



(44) 
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If not for the inhomogeneity in the boundary conditions we could have solved Eas. p5)) . pi)) by convolution with 
the Green's function. Let us denote this convolution U{x,s), i.e. 



U{x,s)^ l\G{x,^,s)^^^d^+ f G{x, s) X + 2 — a, s) u{a, 0). (45) 



As shown in Appendix, inhomogeneity in the boundary conditions introduces two extra terms to the right hand 
side of Eq. ps)) analogous to its last term, so 



U{x, s) = U{x, s) + —G{x, 0, s) w(0, 0) + —G{x, L, s) u{L, 0). 
c c 



(46) 



Since r(a;, ^, t) :— C ^[G{x, ^, s)] properties of Laplace transform immediately imply 



C-'[sGix,^,s)]^Ttix,^,t) 

C''[G{X,C,S)P{C,S)]^ f r{x,^,t-T)p{^,T)dT. 

Jo 



Technically speaking, C-'^[s G{x, ^, s)] = Tt{x,(,t) + r{x,^,0) S{t), but we have ignored the delta function since it 
does not contribute for t > 0. Therefore, solution to the initial-boundary problem Eqs.([T]),(I2]) can be written in the 
form 



u{x, t) 



1 r 



/iiu(0, 0) r(a;, 0, t) + h2u{L, 0) T{x, L, t) + 2hzu{a, 0) T{x, a, t) 



(47) 



1 



rt(x,e,t)u(e,o) + r(x,e,t)u(C,o) 



t i-L 



"'0 



V{x,i,t-T)v{i.T)didT 



If /li 7^ ±1, /ii + /i2 + 2/13 ^ and all roots of algebraic equation Ea. ((T6)) are simple then all our assumptions for 
the validity of Eq. ((32|) are satisfied and we can substitute it into Eq. (|47|) . Taking into account that 



p 00 ^ 



r(x,o,i) = 



hi+h2 + 2/13 



p 00 (k) 

— — s 



fc=ln=-oo 



k=l n= — oo 



A. 
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we have for the vibratory response: 



u(x, t) 



1 



hi + h2 + 2/i3 



hi u{0, 0) + h2 u{L, 0) + 2/i3 u(a, 0) + ^ / u{C, 0) + / r) dr 



E E 

A;— 1 n— — oo 



c/ii u(0, 0) + c/i2 0) ipa{L, s(f )) + 2c ha u{a, 0) (/7a(a, s^f^)] (48) 



Jo 



Formula Eq. (j47p remains valid of course even if the eigenmode expansion is inapplicable. For example, in those 
critical cases when r(x,^,t) is found in closed form, see Eas. pO)) . p^ . we can use it to find the vibratory response 
as well. When substituting into Eg . (ITT)) one should keep in mind that Ht{t — a) = 5{t — a), where 5{t — a) is the 
delta function of Dirac. Its convolution with any function simply returns the function's value at a. Also, when 
terms like — a)H{ct — \x — ^|) are integrated over from to L one has to consider cases x > ^ and a; < ^ 
separately to remove the absolute value, e.g. 



H{i~a)H{ct^\x-^\)w{Od^ 



max{a,a: — ci} J max{a,a;} 



w{0 



Note that the first integral is assumed to vanish if its lower limit exceeds its upper limit, e.g. a > x. Explicit 
expressions are cumbersome and are given in the Appendix. 



9 Numerical issues 

In this section we illustrate the behavior of the system for various values of parameters hi, /i2 and h^. The 
critical cases are of particular interest as they demonstrate somewhat counter-intuitive behavior of the bar. Our 
analytical expressions are calculated using Maple and in some cases we compare results to a MATLAB finite element 
implementation. We subject the bar to initial displacement only in order to illustrate how each damper at the 
boundary and along the bar modifies traveling waves. We use a Gaussian function 0.1 exp(— (x— /i)^/(2 a^))/{a V27r) 
as an approximation of an impulse function where a = 0.1. We first begin with some special values of the parameters 
hi. For all figures c = 1.5 and L — 1.8. 
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9.1 Response of the system for /i2 = 1, = 

Figure [3] shows the response of the system when hi = 0.5, h2 = I and /13 = 0. The right boundary is transparent and 
the left-half of the traveling wave disappears on the right while the right-half get reflected from the left boundary 
and then it too disappears on the right. This behavior can be termed super-stable because the bar comes to rest in 
finite time. Such a phenomenon can only be appropriately observed using explicit solutions and it is the consequence 
of the continuum nature of the system. Once we discretize the system this behavior disappears and the bar could 
only come to rest at infinite time. The reason for this is an exponential nature of the solution for a discrete system. 




Figure 3: Vibratory response for /13 = 0, /ii = 0.5, = 1 (left), and hi = /i2 = 1 (right). 



9.2 Response of the system for hi = h2 = 1, ^ 

Figure 2] depicts the response of the system for hi = /i2 — 1, ^13 = .5. The internal damper is placed at distance 
a — 0.5 from the left end of the bar. The initial pulse consists of two Gaussian functions centered at 0.25 L and 
0.75 L . We see that a wave gets partially reflected by the internal damper. This can be best observed in the Figure 
2 for the left-half of the right impulse. For greater values of the ft,3parameter the internal damper acts almost as if 
it is the fixed point of the bar and waves get reflected to a greater extent while only a small portion of the wave is 
passes through. 



23 



x(m) 



Figure 4: Vibratory response for hi = h2 = 1 and = 0.5 
9.3 Response of the system for /12 = 0.99, /is 7^ 

Figure [S] depicts the response of the system for hi ~ 0.3, /i2 = 0.99 and h^ — 0.7 using 40 eigenfunctions and a 
Gaussian impulse at /i = 0.25 L . The internal damper is placed at distance a — 0.5 L from the left end of the bar. 
In this case the right boundary is almost transparent which precludes formation of standing waves on the right of 
the internal damper. The left boundary, however, is reflective and the standing waves are clearly visible on Figure 
[SI This solution was validated with a finite element implementation (FEM) calculated with 160 elements and the 
error was found to be within 0.001. As the number of elements in FEM was increased this difference became smaller 
indicating the efficiency of the analytical expression. 
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Figure 5: Vibratory response for hi — 0.3, /i2 = 0.99 and — 0.7. 
9.4 Response of the system for hih2 = 1 and = —hi^dis. 

Figure |6] shows the response of the system for hi — 3/10, /12 = 10/3, = —109/60 and a = L/2 with a Gaussian 
impulse at /i = 0.25 L. This is the case 3) from section[5] In this case the system has a double pole at zero and all 
the eigenvalues are imaginary, i.e. there is no damping present in the bar. The amount of energy lost at the left 
and right boundaries is returned into the bar by the damper in the middle. Therefore, the displacement of every 
point on the bar undergoes periodic motion as depicted in Figure [B] Furthermore, this case is significant since FEM 
does not yield the correct result. Our simulations indicate that the internal damper greatly decreases the accuracy 
of FEM eigenvalue computation. As a result, when all the poles are on the imaginary axis errors in the real part 
significantly distort the response. Similar situation occurs in case 4) of section [5] However, in both cases FEM 
response is even more distorted by spurious eigenvalues which we discuss next. 
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Figure 6: u{x,t) for hi = 3/10, /12 = 10/3, /ig = -109/60 and a = 0.5L(left) and (right). 
9.5 Spurious eigenvalues in FEM 

We have encountered a phenomenon for which we do not have an explanation at this time. We have observed 
that for some values of parameters hi the stable continuous system becomes unstable when discretized by the FEM 
method. One set of parameters that produces such a behavior is hi — 0.7, — —1.5 and h^ ~ for which the 
distribution of its eigenvalues is depicted in Figure [71 It is clear that the continuous system is stable since there are 
no eigenvalues with positive real parts. 
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Figure 7: Eigenvalues in the complex plane for hi = 0.7, /i2 



1.5. 



If, however, one discretizes the system an eigenvalue with a positive real part will arise. It can be observed 
that a spurious eigenvalue lies on the positive part of the real axis in the complex plane which forces the discrete 
system to become unstable. The situation is not improved by increasing the number of finite elements. On the 
contrary, the real part of the spurious eigenvalue will increase to infinity making the system even more unstable. 
This is counter-intuitive since we know that the continuous system is equivalent to the discrete one as the number 
of elements tends to infinity. Note that this phenomenon has nothing to do with the system being unconstrained. 
Therefore, it can be concluded that at least some stability regions for some parameters of the continuous system 
become so distorted that after discretization we observe unstable behavior. This phenomenon may have its roots 
in the non-self-adjointness of the continuous system and at this time it is unclear to us how a discretization of such 
a system changes its behavioral pattern. 

10 Conclusions 

We studied longitudinal vibrations of a bar with viscous ends and internal damper. The corresponding eigenvalue 
problem contains the spectral variable in the boundary conditions and has complex- valued, non-orthogonal eigen- 
modes. Behavior of the system is controlled by four dimensionless parameters, the three damping coefficients hi, 
and the relative position of the internal damper a/ L. Our main observations are summarized below. 

Despite the unconventional nature of the eigenvalue problem the eigenmodes can be found explicitly if the 
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eigenvalues are known. When there is no internal damper or it is located in the middle of the bar the eigenvalues 
can be found analytically as well. Otherwise, for rational ajL their determination reduces to solving an algebraic 
equation. Distribution of eigenvalues is hypersensitive to the value of ajL: its complexity grows swiftly with the 
denominator of a/L, and it becomes pseudo-random for irrational ajL. 

When the values of hi are not restricted to be non-negative, i.e. the dampers are allowed to be boosters, effective 
zero damping may occur. Combinations of parameters that lead to zero damping are non-trivial, but can be found 
analytically at least when there is no internal damper or it is located in the middle of the bar. 

Generically, the eigenmodes are sufficient to expand the Green's function. An explicit series solution can then 
be found for the vibratory response, it is a superposition of exponentially damped (or boosted) standing waves. 
These waves however, are complex- valued. Non-generic behavior occurs in two different situations. The first one is 
that the eigenvalues may have higher multiplicity and the associated modes are required for complete expansion. 
This only occurs when Ea. (jl6p has multiple roots or /ii -I- /12 + 2/13 = 0. In practice, one can sidestep the issue by 
slightly perturbing the damping coefficients to resolve the multiplicity. 

The second situation is critical behavior, when the eigenvalues disappear partly or fully. This only happens for 
hi = ±1, but perturbing the coefficients will lead to a qualitatively different picture at large times. When hi = 1 
the corresponding damper turns into a perfect absorber draining all the energy from the bar in finite time. With 
no reflection the standing waves either can not form at all (super-stability) , or can only form on a part of the bar 
accounting for scarcity of the eigenmodes. When hi ^ —1 one of the dampers turns into an infinite amplifier and 
the solution can only exist before the traveling waves reach that damper (super-instability). Mathematically, the 
eigenmode expansion is impossible because the Laplace transformed Green's function is unbounded at infinity and 
is not equal to the sum of partial fractions over its poles. When no eigenmodes exist the solution can be found in 
closed from, but partial cases like /ii ^ 1, /i2 = 1, 7^ require further work. Neither eigenmode expansion, nor 
closed form solution are possible for such cases. This can only happen when the internal damper is present. 

FEM provides good approximation of the vibratory response for small times even in critical cases. However, it 
is unreliable for determination of eigenvalues even when the number of elements is large. Moreover, it universally 
produces eigenvalues with large real parts that do not converge to any actual eigenvalues as the number of elements 
is increased. When these spurious eigenvalues have positive real parts the FEM vibratory response is also unreliable 
for large times. 
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Appendix A: Coefficients of partial fractions 



In this appendix we derive an explicit formula for the coefficients of the partial function expansion of the Green's 
ction. For simplicity, within this ^ 
To this end consider the integral 



function. For simplicity, within this Appendix we write s„, A„ instead of s,^\ Ai*^"* 



(A.1) 



Taking into account that G(^, x, s) and </3a(Cj ^n) both satisfy the boundary conditions the integral / becomes after 
integration by parts, 



/ = (/?a(a;, Sn) + — [hiG{x, 0, s)ipa{0, Sn) + h2G{x, L, s)tpa{L, Sn)] ■ 

c 



(A.2) 



Now we wish to take the limit s — >■ s„. Contrary to intuition, the second term does not tend to when s — )> s„ 
because G has a pole at s„. Indeed, we have from Eq. ((29|) that G{x,S,,s) — + Reg(a;,^, s) where 

Reg(a;,^, s) is regular at s = s„. Substituting this expression we find that 



hi 



+ 0{S - Sn) 



^ ipaix,S„) 



1 + ^ {hiifaiO, SnY + h2^a{L, S„)^) 



(A.3) 



On the other hand, we can compute / in a different way. Since ^'aii^Sn) — 2h3^6{^ — a)(/3a(C, Sn) + ^^aH, Sn) 
from the equation we have 



S S S ~~' S S ' — s 
2h3-6{^ - a)ipa{^, Sn) + -n'Pai^, Sn) " fai^, Sn) = 2/l3 -6{^ - a)lfia{^, Sn) + T-^fa{£., Sn)- 



(A.4) 
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Substituting Ea. (|A.4|) into Eq. (|A.ip yields 



V -^nyS — Sn) I \ ^ n 



S{^-a)ipa{L,Snrd^ + 



A r2 , 



(A.5) 



Combining Eas. (jA.3|) and (jA.5|) finally leads to 



2s, 



c 



h2 



2h. 
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(A.6) 



This formula is derived in 1 with missing boundary terms because the authors overlook that the second term in 
Eq. (jA.2p does not vanish as 5 — y Sn ■ After some algebra and simplifications the formula is given by 



An = --^i 4 (L - a) (f{a, s„)^/i3 + 



2 [L-^aj (hi-l) e-^ + 2 {hi + 1) (l - ] ] ^{a, s„) 



+ i a (/ii - 1)' e-2 ^ - i a {h^ + if e^ ' 



h3 + L{hl-l)}. (A.7) 



Appendix B: Inhomogeneous boundary conditions 



Consider the boundary value problem 



2/13- S{x — a) + — \ U{x, s) ^^2' ^ — ^(^' 



(B.l) 



dU{0,s) hi 

— 1 sC/(0,s) = 71 

ax c 



(B.2) 



dU{L,s) h2 , . 



(B.3) 
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Let U{x, s) be the solution with ji — ^2 — 0. Then U{x, s) — J^^ G{x, ^, s)w(^, s)c?^, where G(x, ^, s) is the Green's 
function of the homogeneous system. We wish to find a solution for general 71, 72. 

The Green's function satisfies L [G(a;, ^, s)] = (5(x — ^) with L — {2h^2^{x — a) + p-^ — and the homogeneous 
boundary conditions. Consider the Green's function of the adjoint problem. Let L* = (^hj,^5{x — a) + — 
be the adjoint operator to L [12, chap 20.3]. The adjoint Green function g satisfies L* [5(2;,^, s)] = 5{x — and 
the adjoint boundary conditions. 

We begin with the integration by parts to obtain 

f-L 

g{x,^,s)L[U{x,s)]dx^[g{x,^,s)U,Ax,s)-g^{x,^,s)U{x,s)]l+ / L*g{x,^,s)U{x,s)dx (B.4) 
^0 

In view of the boundary conditions the first term on the right of Eq. (jB.4[) . the so called surface term, simplifies to 
g{L, ^, 5)72 — g{0, 5)71, while the second term becomes U{£,, s) because of S{x — ^). This yields 

Ui^,s)= [ gix,^,s)L[Uix,s)]dx-g{L,^,s)^2+g{0,^,shi (B.5) 

It is known |12l chap 20.3] that g{x, ^, s) = G(^, x, s). Substituting this into Eq. (jB.5[ ) and interchanging x and ^ we 
finally obtain 

U{x, s) = U{x, s) + jiG{x, 0, s) - 72G(x, L, s). (B.6) 



Appendix C: Vibratory response for hi = h2 = 1 

When hi — h2 — 1 it is possible to obtain a closed form solution for the vibratory response by substituting Eq. p^ 
into Eq. (H71) . Here we give the final expressions for the integrals from the second line of Eq. P?]) . To shorten the 
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formulas we set /(^) := it(^, 0) and g{^) := u{^,0): 



c Jo 



-f{2a - x-tc)H (t- 



2(1 + /is)' 



a — x\ 1 



1 X / a — a; 

f{x + tc)H t 



H{a-x) 



H{x — a) ; 



T{x,i,t)g{i)di 



2(1 + /is) 



/ 5(2a 
Jo 



( a — x 

T 



X — TC)H j 2 ' ~ TC)dT 



2 Jo 



+ ^i^(- + -)^(^--)^-+2(l + 



L_j[^(. + rc)il(r-^Mr 



iJ(a - x) 



+ 



/is 



2(1 + /is) 



y (/(2a — x + Tc)H — — \ J ~^ TCjdr 



\j^9{x- Tc)H i ^-^ - T 1 dr + 



2(1 + 



— /' 

+ /is) io 



g{x — Tc)H ( r ) dr 



H{x — a) ; 



t rL 



Jo 



T{x,^,t-T)p{^,T)d^dT-. 



2c(l + /is) 



rp{^,T)d^d7 

Jo 



+ (1 + /IS)/ ° r p{^,T)d^dT + il + hs) [ ° r p{^,r)d^dT + [ ° [ p{^,T)d^dT 
Jo Jo Jo Jx Jo Ja 



H{a - x) 



+ 



2c(l + /is) 



/ p{^,T)d^dT-hs [ ^ [ pi^,T)d^dr 

Jo Jo Jx 

+ {l + hs) [ " r p{^,T)d^dT + [ " [ p{^,T)d^d7 

Jo Ja Jo Ja 



H{x - a) . 
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